Skip to content

Improve handling of large values by inverse hyperbolic and trigonometric functions#2928

Open
antonwolfy wants to merge 12 commits into
masterfrom
update-inverse-trig-functions
Open

Improve handling of large values by inverse hyperbolic and trigonometric functions#2928
antonwolfy wants to merge 12 commits into
masterfrom
update-inverse-trig-functions

Conversation

@antonwolfy
Copy link
Copy Markdown
Contributor

@antonwolfy antonwolfy commented May 20, 2026

The dpnp.tensor.acosh function was returning incorrect results (infinity) for complex numbers with large negative real parts (see #2924).

The issue was caused by the threshold used 1/epsilon to decide when to switch from acos in SYCL's experimental complex extension (exprm_ns::acos()) to a log-based formula.
The SYCL exprm_ns::acos() returns incorrect results for large negative real parts, due to used formula z + sqrt(z² - 1), which leads to log(0) = -infinity.

The PR proposes to changed the threshold from 1/epsilon to sqrt(1/epsilon)/2 for all inverse trig/hyperbolic functions. This is the precision loss point where sqrt(z² ± 1) calculations become unstable.

Also the PR reduces code duplication by creating shared casinh(), casin(), catanh(), catan() functions in new complex_math.hpp header.

  • Have you provided a meaningful PR description?
  • Have you added a test, reproducer or referred to an issue with a reproducer?
  • Have you tested your changes locally for CPU and GPU devices?
  • Have you made sure that new changes do not introduce compiler warnings?
  • Have you checked performance impact of proposed changes?
  • Have you added documentation for your changes, if necessary?
  • Have you added your changes to the changelog?

@antonwolfy antonwolfy added this to the 0.21.0 release milestone May 20, 2026
@antonwolfy antonwolfy self-assigned this May 20, 2026
@github-actions
Copy link
Copy Markdown
Contributor

View rendered docs @ https://intelpython.github.io/dpnp/pull/2928/index.html

@antonwolfy antonwolfy linked an issue May 20, 2026 that may be closed by this pull request
@coveralls
Copy link
Copy Markdown
Collaborator

coveralls commented May 20, 2026

Coverage Status

coverage: 78.239%. remained the same — update-inverse-trig-functions into master

@github-actions
Copy link
Copy Markdown
Contributor

github-actions Bot commented May 20, 2026

Array API standard conformance tests for dpnp=0.21.0dev0=py313h509198e_54 ran successfully.
Passed: 1356
Failed: 4
Skipped: 16

@antonwolfy antonwolfy marked this pull request as ready for review May 20, 2026 15:55
@ndgrigorian
Copy link
Copy Markdown
Collaborator

The PR proposes to changed the threshold from 1/epsilon to sqrt(1/epsilon)/2 for all inverse trig/hyperbolic functions. This is the precision loss point where sqrt(z² ± 1) calculations become unstable.

should this be filed as an issue with DPC++ compiler? Since most implementations should have logic to prevent overflow/infs in such cases

@antonwolfy
Copy link
Copy Markdown
Contributor Author

should this be filed as an issue with DPC++ compiler? Since most implementations should have logic to prevent overflow/infs in such cases

Yes, we probably should, but I wonder that the extension still in the experimental state.
And we already had such branch previously with 1/epsilon threshold, probably due to the exact the same reason.

I also tried to implement the kernel through SYCL's C complex math functions since marked as supported extension:

// Forward declare SYCL's C complex math functions with C linkage
extern "C" {
extern __DPCPP_SYCL_EXTERNAL_LIBC float __complex__ cacosf(float __complex__ z);
extern __DPCPP_SYCL_EXTERNAL_LIBC double __complex__ cacos(double __complex__ z);
}
...
        if constexpr (is_complex<argT>::value) {
            using value_type = typename argT::value_type;

            if constexpr (std::is_same_v<value_type, float>) {
                auto result = ::cacosf(*reinterpret_cast<const float __complex__*>(&in));
                return *reinterpret_cast<const std::complex<float>*>(&result);
            }
            else {
                auto result = ::cacos(*reinterpret_cast<const double __complex__*>(&in));
                return *reinterpret_cast<const std::complex<double>*>(&result);
            }
        }
...

but faced exactly the same #2924 issue, since the SYCL does provide cacos implementation, but with exactly the naive formula:

double __complex__ w = clog(z + csqrt(__sqr(z) - 1.0));

Btw I'd propose to merge that PR, in background I'll submit the issue towards the compiler.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

acosh returns wrong result

3 participants